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We examine the origins of azimuthal correlations observed in high energy proton-nuclens colli¬ 
sions by considering the simple example of the scattering of uncorrelated partons off color fields in a 
large nucleus. We demonstrate how the physics of fluctuating color fields in the color glass conden¬ 
sate (CGC) effective theory generates these azimuthal multiparticle correlations and compute the 
corresponding Fourier coefficients within different CGC approximation schemes. We discuss in 
detail the qualitative and quantitative differences between the different schemes. We will show how 
a recently introduced color field domain model that captures key features of the observed azimuthal 
correlations can be understood in the CGC effective theory as a model of non-Gaussian correlations 
in the target nucleus. 
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I. INTRODUCTION 

Azimuthal anisotropies of multiparticle correlations 
observed in small systems such as those produced in 
p-|-Pb, d-|-Au, or ^He-|-Au collisions have been computed 
in various theoretical frameworks. Within different cal¬ 
culations these correlations are either dominantly due to 
initial state parton correlations in the projectile and tar¬ 
get [1-15] or from final state correlations that are gen¬ 
erated by the collective flow of matter produced in the 
collision [16-22]. Since both of these frameworks are able 
to describe key features of the data, disentangling the 
different effects and understanding how correlations are 
generated in both approaches is essential to obtain novel 
insight into the QCD dynamics of ultradense parton sys¬ 
tems. 

We will focus in this work on multiparticle correla¬ 
tions which are generated in the initial state. We will 
discuss within the color glass condensate (CGC) effec¬ 
tive theory of high energy QCD [23] the origin of these 
correlations and we will critically examine the assump¬ 
tions underlying different calculations in this framework. 
Even though some of the models may appear very dif¬ 
ferent, their common features and their differences can 
be understood systematically as approximations to the 
underlying QCD dynamics. 

We will orient our discussion of initial state correlations 
within the “dilute-dense” power counting in the CGC, 
where the incoming parton densities are assumed to be 
small in the projectile but large in the target nucleus. 
In this limit, analytical and numerical computations 
are comparatively simple and thus permit systematic 
comparisons between different approximation schemes. 
We note however that the kinematic region where the 
strongest azimuthal correlations are seen in experiments 
correspond more to a “dense-dense” situation as the par- 
ton densities are also large in the incoming projectile. A 
systematic power counting then requires one to solve clas¬ 
sical Yang-Mills equations in the presence of projectile 


and target color sources, which has only been achieved 
numerically [14, 24]. While there are important quali¬ 
tative differences between dense-dense and dilute-dense 
systems, we believe that the lessons inferred from ana¬ 
lytical and numerical studies the dilute-dense case can 
nevertheless be valuable for the discussion of the phe¬ 
nomenologically more relevant dense-dense collision sys¬ 
tems. 

We will here compute the two particle correlation func¬ 
tion for quarks scattering off a large nucleus in the 
dilute-dense CGC description. The main ingredient of 
this computation is the so-called “dipole-dipole” corre¬ 
lator of light-like Wilson lines and we will compute this 
quantity in the following approximations: i) the Gaus¬ 
sian two gluon exchange (“Glasma graph”) approxima¬ 
tion [1-5, 25-29] and ii) the nonlinear Gaussian approx¬ 
imation [30-40]. We will determine the second, third, 
and fourth azimuthal Fourier coefficients of the two par¬ 
ticle correlation function within these two approxima¬ 
tion schemes and compare the results to numerical lat¬ 
tice simulations of the full correlator in the McLerran- 
Venugopalan (MV) model [41-43] and after renormal¬ 
ization group (JIMWLK) evolution [44-47] of the MV 
model initial conditions to higher rapidities [13]. We 
will further discuss how our computations relate to the 
“color field domain model” introduced in [10, 11, 48, 49] 
based on ideas developed in [6, 7]. Since many features 
of this model appear similar to those discussed previ¬ 
ously [1, 2, 4, 5], it is important to understand the in¬ 
terpretation and justification for this model from first 
principles. We will demonstrate that the effects of the 
color field domain model can be reproduced in the CGC 
effective theory if non-Gaussian correlations are assumed 
to play an important role. 

This paper is organized as follows. In the next sec¬ 
tion, we shall discuss the physical picture of how initial 
state multiparticle correlations are generated and derive 
the leading order expressions for single and double inclu¬ 
sive distributions. These can be expressed in terms of 
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FIG. 1: (Color online) Color electric fields inside the nucleus 
fiuctuate on an event by event basis. 


correlators of lightlike Wilson lines, which we will calcu¬ 
late in Sec. 3 within the Glasma graph and nonlinear 
Gaussian approximations. We then compute the Fourier 
moments Vn of two particle correlations in both approx¬ 
imation schemes in Sec. 4 and compare our results with 
numerical lattice simulations of the full correlation func¬ 
tion. In Section 5 the analytical and numerical results 
obtained are then compared to the color field domain 
model. We first cast our analytical results in terms of 
color electric field correlators and make a direct com¬ 
parison with expressions for the same correlator in the 
color field domain model. We will demonstrate that our 
results provide a clear interpretation of the color field 
domain model which clarifies the discussion in the recent 
literature. We end with a summary of the results of the 
paper and an outlook on further research directions in 
computations of multiparton correlations in high energy 
QGD. 


II. MULTIPARTICLE CORRELATIONS EROM 
ELUCTUATING COLOR FIELDS 


We begin our discussion of the physics of initial 
state correlations with the simplest possible example of 
the high energy scattering of individual (uncorrelated) 
quarks off a large nucleus. Our general picture is that 
each parton scatters independently off the color field of 
the nucleus receiving a transverse momentum kick in the 
process. As noted previously [6, 7, 25], the color fields 
fiuctuate from event to event and are locally organized 
in domains of size ~ 1/Qs as illustrated in Fig. 1. When 
two (or more) quarks scatter off the same domain, they 
will receive a similar kick whenever they are in the same 
color state. This leads to a correlation which is sup¬ 
pressed by 1/Ac^ (in the limit of large N^) and the num¬ 
ber of domains QsS±, where S± denotes the transverse 
area probed by the projectile. We will now discuss this 
physical picture in in more detail and further develop its 
quantitative implementation along the lines of the dis¬ 
cussion in Ref. [13]. 


A. Single quark scattering 


Within the CGG formalism, the color fields inside the 
target nucleus are determined by the solution of the clas¬ 
sical Yang-Mills equations 

[77^,Fn=J", (1) 

where the eikonal current is given in terms of the 
density of color charges p inside the target nucleus as 

= 5^~p{-x.,x'^) . ( 2 ) 


The solution to the classical Yang-Mills equations takes 
the well known form [50] 


A (x,x~^) 


P(x,a:+) 


( 3 ) 


where = didi is the 2-dimensional Laplacian. The 
scattering of an incoming quark inside the projectile can 
be described to leading order accuracy in by the so¬ 
lution of the Dirac equation 

(i]/> — m)'!' = 0 , (4) 


in the presence of the background field of the target in 
Eq. (3). One finds that the forward scattering amplitude 
of a quark with momentum p to scatter off the color fields 
in the target is given by ^ 


(out,q|in,p) 


where 


V{x) — V exp 



( 5 ) 

( 6 ) 


denotes the Wilson line at a spatial position x in the 
fundamental representation. 

Within the leading order dilute-dense framework, it 
is then straightforward to compute the single inclusive 
distributions of quarks in a high energy projectile after 
scattering from a nuclear target 


dNg 

d^p 


(out,p|/5|out,p) , 


( 7 ) 


where p is the reduced one particle density matrix in the 
probe. This general expression can be rewritten explicitly 


^ Since we are primarily interested in the transverse coordinate 
dependence, we have omitted a delta function for longitudinal 
momentum conservation as well as the spin structure to lighten 
the notation. We refer to [51] for the complete expression. 

^ Our expression generalizes the one given in [51] by replacing 
the collinear quark distribution with a Wigner function Wg(b, k) 
that is a function of both the k of quarks in the projectile and 
their impact parameter b. Equivalent expression for gluons, dif¬ 
fering only by the representation of the Wilson lines, have ex¬ 
plicitly been derived in [31, 52]. 





3 


dNg 

d^p 



(2^ 


W^,(b,k)<p(b,p 


k). (8) 


The Wigner function Wq{h,k) characterizes the trans¬ 
verse momentum and position distribution of incoming 
quarks inside the projectile and is defined to be 


Wq{h,k) 


d2q 


in, k -I- 


q 

2 



g-iqb 


(9) 

For the illustrative purpose of this paper, it is sufficient 
to choose a Gaussian form 


lF,(b,k) 


— e"'" 

TT^ 


( 10 ) 


with a dimensionful constant B characterizing the trans¬ 
verse area of the projectile. The dynamics of interest to 
us is given by the unintegrated gluon distribution of the 
target nucleus 


<P(b,k)=y d^r £) (b-b ^,b-I) , (11) 

which represents the distribution of momentum transfers 
from the target that contribute to the corresponding mo¬ 
mentum distribution of the scattered quark. Here 

i?(x,y) = (I?(x,y)) (12) 

is the expectation value of the dipole operator 

I?(x,y) = i-Tr [H(x)Ht(y)] , (13) 

which results from the product of the forward scattering 
amplitude in Eq. (5) and its complex conjugate equiv¬ 
alent, to obtain the single inclusive probability [53-55]. 
Written fully in coordinate space our expression for the 
single inclusive multiplicity (8) has the form 


dA^g _ 1 1 

d^p ttB (27r)2 


d'^Kd^yD{x,y) 


gip.(x-y) 


e~^ 



(14) 


which is recognizable as the one used in [13] up to a con¬ 
stant factor which cancels in the anisotropy coefficients 
Vn- 

Since we are interested in the scattering of a small 
probe - such as a quark inside a proton - off a large 
nucleus, it is reasonable to neglect the impact parameter 
b = (x -b y)/2 dependence in the target. Because on 
average there is no preferred direction in the transverse 
plane, the expectation value Z?(x, y) then depends only 
on the magnitude of the transverse separation r = x — y 
of the dipole^. 


D{x,y) = D{\x-y\) . (15) 


® The dipole expectation value can in general depend on the rota- 


Since the dipole coordinate r is the conjugate variable 
to the momentum transfer p — k for a single quark 
scattering, the symmetry in Eq. (15) ensures that the 
momentum transfer to each individual quark is on 
average symmetric with respect to the azimuthal angle. 


B. Double inclusive spectrum and multiparticle 
correlations 

We will now consider the case where two quarks scat¬ 
ter independently off the same nucleus and shall study 
the correlations between the two scattered quarks. We 
will make the simplifying assumption that the momenta 
of the two incoming quarks are initially uncorrelated ^ 
- the two particle distribution VFg,(bi, ki, b 2 , k 2 ) of in¬ 
coming quarks factorizes into the product of single quark 
distributions. 


Hjjq(bi,ki,b2,k2) = lT,(bi,ki)Hjj(b2,k2) . 


(16) 


The double inclusive distribution of scattered quarks 
then takes the form 


d^Af 


d^pi d^p:, 


= / d"bid"b2 


d^ki 


d^ko 


(27r)2 J (27r)2 


X Wq(bi,ki)W,(b2,k2) 

d^ri d^r2e 


i(pi-ki).ri gi(P2-k2)'r2 



xv(h2 + Y>b2 - 



(17) 


While we assumed the transverse momenta ki and k 2 
of the two incoming quarks to be uncorrelated, one 
observes from Eq. (17) that this is no longer the case 
for the momenta pi and p 2 of the scattered quarks. 
Since both quarks scatter off the same nucleus, the 
momentum transfers pi — ki and p 2 — k 2 are correlated 


tional invariants |b|,|r| and b • r. Once one performs the averag¬ 
ing over the b distribution of the projectile, only the dependence 
on r remains. We note however, that such an averaging can 
effectively introduce non-Gaussian correlations between several 
dipoles. These can affect multiparticle correlations at a charac¬ 
teristic momentum scale given by the inverse impact parameter 
dependence. 

We are mostly concentrating on the near-side “ridge” correlation 
for semihard momenta ~ Qs. Thus we are neglecting back-to- 
back momentum correlations [56, 57] that are particularly impor¬ 
tant for the away-side “jet” peak [3, 58, 59] and contributions at 
the nonperturbative small intrinsic transverse momentum scale 
of the probe. Such correlations should be taken into account in 
a full comparison with data. 
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with each other, giving rise to azimuthal correlations of 
the scattered quarks. 

We note that the above model can be generalized in 
a straightforward way to study correlation functions in¬ 
volving more than two particles. Such higher order corre¬ 
lation functions can be related to higher order correlation 
functions of dipole correlators. For example, four quark 
correlations in this model involve expectation values of 
products of four dipoles in the fundamental representa¬ 
tion. 


III. DIPOLE-DIPOLE CORRELATOR 


The discussion in the previous section shows that 
all the features of two-particle correlations are encoded 
in the expectation value of the dipole-dipole correlator 
(I?(x, y)I?(u, v)). We will now study the properties of 
this correlator and discuss approximation schemes that 
have been frequently employed in the literature. 


A. Glasma graph approximation 


The basic properties of the two particle correlation can 
be understood by studying the interaction between the 
incoming quark and the target nucleus in terms of multi¬ 
ple gluon exchanges. Clearly the dynamics of these gluon 
exchanges depends on the nature of color charge corre¬ 
lations in the target and needs to be specified. A simple 
model of such correlations is the MV model [41-43]. In 
this model, the underlying distribution of color charges in 
the nucleus is assumed to be a Gaussian distribution such 
that all multigluon correlations are uniquely determined 
by the two gluon correlation function. 

A further approximation that simplifies the computa¬ 
tion considerably is to assume that each of the quarks 
in the projectile exchanges only two gluons with the tar¬ 
get nucleus. We will refer to this combination of the two 
gluon exchange approximation and Gaussian statistics as 
the Glasma graph approximation, a term first introduced 
in Ref. [2]. This approximation has been used in a num¬ 
ber of phenomenological studies of ridge correlations in 
high energy collisions [1, 2, 4, 5, 25-29]. 

More specifically, the Glasma graph approximation 
can be understood by examining the expression for the 
dipole-dipole correlator in terms of the color field. When 
the density of color charges in the target nucleus gp is 
small - corresponding to a dilute-dilute situation - one 
can perform an expansion of the path ordered Wilson 
line around the identity matrix, representing an expan¬ 
sion in terms of the number of gluons exchanged between 
the projectile and the target. In order to keep the nota¬ 
tion as light as possible, we will denote the Wilson lines 


as F(x) = exp(—iA(x)) in the following®, so that the 
expansion takes the form 

V(x) l-zA(x)-iA2(x) (18) 

+ ^A3(x) + ^A4(x) + ... . 

Evaluating the color traces for A(x) = A“(x)t“ as 

ro& 1 

tr[rt®] = — , tT[eth^] = -f , (19) 

the dipole operator to ©(A®) within this approximation 
takes the form 


I?(x,y)-1~ 

1 

~ 2^ 


i 

64W 


(a“ A® - 2A“ A® + A“ A®) 

(a“ A® A^ - 3A“A® A= 

+ 3A“ A® A= - A“ A® A^) + ... 


• ( 20 ) 


In the Gaussian approximation for the correlations of A 
the C and V odd C>(A®) term vanishes in the expectation 
value of the dipole operator and only the 0{A^) remains. 
We will denote the AA correlator as 


(Aa(x)Ab(y)) = 5 “'’ 7 (x - y) , (21) 

which defines the correlation function 7 (r). With this 
definition, we obtain [53] 

£)(x - y) ~ 1 - Cf ( 7 ( 0 ) - 7 (x - y)) . 

Following the same logic as previously for the dipole 
expectation value, we will now compute the dipole-dipole 
correlator by expanding up to 0{A^) in the coupling con¬ 
stant. The assumption of Gaussian statistics means that 
the four point correlation function can be expressed in 
terms of the two point function in Eq. (21) as 

(Aa(x)At,(y)Ac(u)Ad(v)) = (22) 

jaf>^cd 7 (x _ y)^(u - v) 

+ , 5 “( 5 “ 7 (x - u) 7 (y - v) 

-I- _ v)7(y — u) . 

Thus one obtains 


(I?(x, y)V{u, v)) cs Dix - y)7^(u - v) 

(x — u) — 7 (x — v) 


(-<2 


2(A^c" - 1) 


- 7(y - u)-H 7(y - v)) , (23) 


^ Careful path ordering -will not modify the results obtained as 
long as the correlations are local in x"*". 
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Disconnected contributions 




One finds that the terms in Eq.(24) that do not de¬ 
pend on all four coordinates x, y, u, v only contribute 
for very soft momenta p, k on the order of the incom¬ 
ing momentum of the projectile ~ l/'s/S- Since 

we are focusing on the dominant semihard gluons with 
Pt, Qt ^ Qs ^ ^IVB we will neglect these contributions 
in the following and replace Eq. (24) with 


(V(x, y)V(u, v)) ~ D(x - y)D(u - v) 

^ ^(£'(y-u)i:>(x-v) 


-f 


(Jv/ - 1 ) 


+ D(y-v)D(x-u)y (25) 


Equation (25) is in fact the form used in the Glasma 
graph papers [1, 2, 4, 5, 25-29]. In particular, in the 
comparisons to data performed in [2, 4, 5], the expec¬ 
tation value of the dipole correlator is taken to satisfy 
the Balitsky-Kovchegov (BK) equation [60, 61]. Thus 
the Glasma graph approximation, as employed in phe¬ 
nomenological computations, corresponds to the merging 
of gluon ladders from multiple sources (all localized on 
a transverse scale ~ 1/Qs) into a single ladder exchange 
represented by the dipole correlator. 


FIG. 2: (Color online) Classes of diagrams contribnting to 
the dipole-dipole correlator in double-gluon exchange approx¬ 
imation. 


B. Nonlinear Gaussian approximation 


which, to this order of the approximation, is consistent 
with 


(V(x, y)V(u, v)) ~ D(x - y)D(u - v) 


1 


2(A^c" 


— (d(x-v)-£I(x-u) 
-i:i(y-v)-b£i(y-u)) . 


(24) 


One observes from Eq. (24) that the leading Nc contri¬ 
bution to the dipole-dipole correlator factorizes into the 
product of single inclusive averages - corresponding to 
the independent scattering of two quarks. Diagrammati- 
cally this corresponds to the disconnected contribution in 
Fig. 2. Genuine correlations are contained in the second 
term of Eq. (24) and suppressed by a factor of 1/ (— 1) 
as pointed out previously in [25, 27, 28]. These corre¬ 
spond diagrammatically to the connected graphs in Fig. 2 
and feature all possible contractions between the x, y and 
u, V dipoles. 

When calculating two particle correlations using 
Eq. (17), one needs a Fourier transform of the dipole- 
dipole correlator in Eq. (24) to obtain the momentum 
transfer p — k from the interaction with the target. 


If we restrict ourselves to Gaussian correlations of 
gluon fields in the target, it is possible to resum the 
multiple gluon exchange contributions to the dipole am¬ 
plitude and the dipole-dipole correlator analytically to 
all orders®. Generalizing the notation from the glasma 
graph case, one obtains [53] 


D(x - y) = exp (^Cp (jix - y) - 7 ( 0 ))) . (26) 


This expression is a part of the nonlinear Gaussian ap¬ 
proximation because it includes all orders in A, evaluated 
with a the Gaussian AA correlator in Eq. (21). This is 
to be contrasted with the two-gluon exchange approxi¬ 
mation Eq. (22) which was expanded to the lowest order. 
Using a well known algorithm [31, 32, 34-36, 38, 38- 
40, 55] for computing higher point correlators one can 
obtain the dipole-dipole operator expectation value to 


As suggested by our previous discussion, these become multiple 
ladder exchanges upon BK or JIMWLK evolution of the dipole 
and dipole-dipole correlators. 
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all orders in A assuming a Gaussian AA correlator as [62] 


(I?(x, y)I?(u, v)) = i^(x - y)D{u - v) 

^F(x,u;y,v) + TA i^(x,y;u,v) 


^^/A 


2y/A 


iVc^VA 


/ f(x,u;y,v) - y/A _ f(x,y;u,v) \ 

2yA N.'^VA ) 



where -F(x, y; u, v) is defined to be 


X 

u 



y -X' 

V u 




y 

1/ 


FIG. 3: (Color online) Example of a higher order contri¬ 
bution to the non-linear Gaussian approximation, which con¬ 
tributes to the X y anti-symmetric part of the dipole-dipole 
correlator at C)(A®). 


, . 1 , f D(x — u.)D{y — \-)\ , . 

F(x,y;u,v) = —In —-—--- . (28 

Cf \D{x-v)D{y-u)J 

Here A is short for A(x, y; u, v) and is given by 

A(x, y; u, v) = (x, u; y, v) 

+ y; v)F(x, v; u, y) . (29) 


We will refer to Eq. (27) as the “nonlinear Gaussian ap¬ 
proximation” for the dipole-dipole correlator. It is exact 
in the McLerran-Venugopalan (MV) model where dipole 
and multipole correlators depend nonlinearly on the A 
fields, but the correlators of A’s are Gaussian. One can 
easily check that the two-gluon exchange limit of the non¬ 
linear Gaussian (27) is the same as the two-gluon ex¬ 
change approximation (23). 

One can obtain sight into this relation by taking the 
large limit for a constant D{x — y): 


operator in Eq. (20). Note that while the expectation 
value of the odd term is zero in the Gaussian approxi¬ 
mation, the expectation value of its square is not, but is 
proportional to d^^bc^abc ^ xhe 

contributions from the odd terms correspond diagram- 
matically to the processes depicted in Eig. 3 and take the 
form 


(F(x,y)I?(u,v)) - (I?(y,x)I?(u,v)) ~ 


Cl Vc^-4 

(Ae^ - 1)2 


X - u) - 7(y - u) 


-7(x 


v)+7(y-v)) . 


(31) 


While for Nc = 2 the antisymmetric contribution van¬ 
ishes identically to all orders, it is nonvanishing for 
Nc > 3 and suppressed by a factor of 1/Ac^ relative 
to the disconnected contribution in the large Nc limit. 


(I?(x, y)V{u, v)) = D{x - y)D{u - v) 

2 


1 


In 


In 


£>(x-u)n(y-v) ' 
^ n(x-v)n(L 

'£>(x-y)£>(L 


D{x-v)D{u-y) 


+D{x—y)D{u—'v) In 


n(x-v)n(u-y) 

£>(x-y)Fi(u-v) 


£)(x-v)C(u-y) 


- 1 


(30) 


which shows again that the leading Nc contribution cor¬ 
responds to the independent scattering of two quarks, 
while genuine correlations are Nc suppressed. 

While the nonlinear Gaussian approximation repro¬ 
duces the double gluon exchange (Glasma graph) ap¬ 
proximation at C(A^) in the dilute limit, it also con¬ 
tains a series of higher order terms. An important sub¬ 
set of higher order contributions corresponds to the di¬ 
agrams which separately break the x — y —> y — x and 
u — V —V — u symmetries. As we will discuss shortly 
these contributions are responsible for generating the odd 
moments (ua, us,...) in the Fourier expansion of the corre¬ 
lation function. One finds that the leading contribution 
to the X ^ y antisymmetric part can be associated with 
the square of the C and V odd contribution to the dipole 


IV. AZIMUTHAL CORRELATIONS IN QUARK 
NUCLEUS SCATTERING 

We will now discuss the azimuthal correlations of 
quarks scattering off a large nucleus and present com¬ 
parisons of results for these correlations within different 
approximation schemes introduced in the previous sec¬ 
tions. To further quantify the correlations introduced 
by the scattering, we will decompose the double inclu¬ 
sive distribution in Eq. (17) into Fourier modes in the 
relative azimuthal angle A(j) between the two scattered 
quarks, 

d.^ 

—-vy—oc 1-b V2V„A(Pi,P2)cos(nA(;i). (32) 

d^pi d^p2 

The familiar coefficients u„{2}(pt) can be obtained from 
the V„A using [63] 

roll' X G„a(pt,Pt®0 
Vn{2}{pT) = , =, (33) 

yl4A(PT®^^’T''0 

where represents a range in transverse momentum 
corresponding to an experimental reference bin. 
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A. Analytic estimates 


B. Numerical results 


Before we turn to a discussion of numerical results, it 
is useful to obtain further analytic insight into the cor¬ 
relation functions themselves by considering the limit of 
a large probe with small intrinsic transverse momentum 
{QgB 1). With the Gaussian Wigner distribution in¬ 
troduced in Sec. II it is convenient to absorb a part of 
the Wigner distribution into the definition of a modified 
dipole distribution 


7(P) = 


d^r 


i D{r) 


(34) 


such that the single inclusive distribution ( 8 ) becomes 


dNg 

d^p 


I 


7(P) • 


(35) 


Within the Glasma graph approximation the double in¬ 
clusive distribution can then be evaluated by combining 
Eqs. (17) and (25) as 


1 


(PN 

d^pd^q (271)"* 
1 


7(p)7(q) 


(A^c" - 1) 


g-(p+q)"B/2 ( ^ 1 P ^ 




which further reduces to 


p + q 


(36) 


r7(p)7(q) 


pn 

d^pd^q (271)"* 

2^ (j(2)(p + q)+j(2)(p_q)) 


1-f 


B{N,^ - 1) 


(37) 


in the limit of a large probe or at high momenta, where 
the intrinsic transverse momentum of the probe can 

B 1,2 

be neglected and we can approximate 2TrBe~^‘^ —>■ 

(27r)^(J^^)(k). While the first term inside the square 
bracket corresponds to the disconnected contribution 
and does not contain any correlations, the connected 
term gives rise to a correlation which is suppressed by 
1/(W" - 1) and l/{QlB). 

One also observes from Eqs. (36) and (37) that the two 
particle correlation function has a similar structure as the 
collinear limit of the Glasma graph computation [ 1 , 2, 4, 
5, 25-29]. It features a near side (p « q) contribution 
as well as one on the away side (p « ~q)) which in 
terms of the Fourier decomposition in Eq. (32) give rise 
to even harmonics V 2 ,V 4 ,.... Odd harmonics U 3 , U 5 ,... are 
not present in Eqs. (36) and (37) as the above expressions 
are manifestly symmetric under p — > — p. 


In order to establish the quality of the different ap¬ 
proximations, we will now compare the results for the 
Fourier coefficients to numerical lattice computations 
that fully evaluate the correlation functions of Wilson 
lines. We begin with a comparison in the McLerran- 
Venugopalan (MV) model where the Wilson lines are 
generated from a Gaussian ensemble of fluctuating color 
charges. Following the numerical procedure of Ref. [64], 
a set of Wilson line configurations is generated according 
to Eq. ( 6 ); these are then employed to extract numeri¬ 
cally the expectation value of the dipole operator D{r). 
With the expectation value of the dipole operator we can 
then compute the double inclusive spectrum in Eq. (17) 
using the dipole-dipole correlator in the Glasma graph 
approximation of Eq. (25) and in the nonlinear Gaussian 
approximation of Eq. (27). We also compute the double 
inclusive spectrum directly from the lattice Wilson lines 
using the procedure described in Ref. [13]. For each of 
these three different double inclusive spectra, we deter¬ 
mine the Fourier coefficients using Eqs. (32) and (33). 
We choose the reference momentum to be = px such 
that Vnipx) = \/VnAiPT). Our results for the Fourier co¬ 
efficients Vnipx) in the MV model are shown in the left 
panel of Fig. 4. 

We also perform, as discussed in [13], the JIMWLK 
rapidity evolution of the Wilson lines for y = 7.6 units in 
rapidity for SU(3) and for y = 12.4 units for SU(2) with 
the same running coupling formula and an initial satu¬ 
ration scale Qs/Aqcd = 3.7 in both cases. We use the 
running coupling prescription for the JIMWLK equation 
proposed in Ref. [65]. We then compute again the dou¬ 
ble inclusive spectrum and the Fourier harmonics directly 
from the lattice Wilson lines, as well as from the lattice 
result for the dipole. The results including JIMWLK 
evolution are shown in the right panel of Fig. 4. For 
the MV model the probe size is BQ^ = 3.7 and for the 
JIMWLK simulations BQ^ = 2.5; for a discussion of the 
R-dependence see [13]. 

We note that for the MV model case the nonlinear 
Gaussian approximation agrees perfectly with the di¬ 
rect numerical calculation for all Vn up to the numeri¬ 
cally accessible values of px- This is of course a sim¬ 
ple numerical check of the analytical expressions and 
the agreement should be exact because non-Gaussianities 
are absent by definition in the MV model. In contrast, 
the Glasma graph approximation deviates significantly 
from the numerical result. As discussed previously, the 
Glasma graph result does not have any odd harmonics. 
However the even Fourier coefficients V 2 and 774 too show 
significant deviations from the exact numerical result es¬ 
pecially around pt ~ SQsj demonstrating that nonlinear¬ 
ities are significant even at larger px. The Glasma graph 
approximation should be exact in the high momentum 
limit when jpj, jqj, Jp —q] and |p + q| are all large. How¬ 
ever, even for a large jpj = jqj the azimuthal harmonics 
receive contributions from |p ± q| < Qs where the two- 
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FIG. 4: (Color online) Azimuthal correlations v„ for scattering of two independent quarks off a large nucleus in the MV model 
(left) and after JIMWLK evolution (right). Solid lines correspond to results from the lattice simulation without additional 
approximation, dash-dotted lines (with squares) show results in the non-linear Gaussian approximation, and dotted lines (with 
circles) correspond to the Glasma graph approximation. 




FIG. 5: (Color online) Nc dependence of the azimuthal correlations v„, scaled by the color factor \/ Nc'^ — 1 for SU(2) and 
SU(3) gauge theory. The results are computed using the numerical lattice calculation for the MV model (left) and after 
JIMWLK rapidity evolution (right). 


gluon exchange approximation is not very accurate. 

With JIMWLK rapidity evolution, the distribution of 
color charges is no longer explicitly Gaussian. There¬ 
fore, although non-Gaussian contributions were not seen 
in the operators studied previously in [37], it is possible 
that JIMWLK evolution of the azimuthal anisotropies 
will introduce non-Gaussian contributions. Indeed this 
is seen in Fig. 4 (right) where we observe a deviation of 
the nonlinear Gaussian approximation from the numeri¬ 
cal result from solving the JIMWLK equations. Further¬ 
more we see that both the numerical and the nonlinear 
Gaussian results for all are reduced by the JIMWLK 
evolution. On the contrary, the Glasma graph results, 
after JIMWLK rapidity evolution of the Wilson lines, 
are roughly the same as those for the MV model. Thus 
while better agreement of the Glasma graph approxima¬ 
tion with the nonlinear Gaussian and full JIMWLK re¬ 
sults is seen after rapidity evolution, this agreement is a 


fortuitous numerical coincidence. 

We have also analyzed the W dependence of our re¬ 
sults for Vn employing the full numerical calculation of 
the dipole-dipole correlator. Our results for the SU(2) 
and SU(3) gauge theory are shown in Fig. 5, where we 
scale the Fourier coefficients V 2 and V 4 by the color fac¬ 
tor y/— 1. This is because azimuthal correlations in 
the double inclusive spectrum contain an overall factor 
of l/(iVc^ — 1) (see e.g. Eq. (37)) and the u„’s are related 
via a square root to the Fourier coefficients VnA in the 
expansion of the double inclusive spectrum. In fact, we 
find that this scaling works nearly perfectly both in the 
MV model case (Fig. 5 (left)) and the JIMWLK evolved 
case (Fig. 5 (right)); in the former, small differences are 
seen in V 2 and V 4 for only for large pt, where lattice cutoff 
effects can already have an effect. 

Finally, we demonstrate the dependence of our results 
on different choices for the reference transverse momen- 
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FIG. 6: (Color online) Comparison of azimuthal correlations V 2 in three different approximations for different choices of the 
reference momentum in the MV model (left) and after JIMWLK evolution (right). We compare pr = to the case of a 
wider reference momentum bin 0.5 Qs < Pt‘^^ < 3Qa- Solid lines correspond to results from the lattice simulation without 
additional approximation, dash-dotted lines (with squares) show results in the non-linear Gaussian approximation, and dotted 
line (with circles) correspond to the Glasma graph approximation. 


turn in Fig. 6. This is an interesting exercise because 
similar studies can be performed with the experimen¬ 
tal data and will help to distinguish between different 
models. We find a clear suppression of the signal to the 
previously regarded case = pp when employing a 
fixed reference bin 0.5 Qs < < ^Qs- One observes 

that for the Gaussian correlations in the MV model, this 
effect is particularly strong at large px in the Glasma 
graph approximation. In case of the JIMWLK-evolved 
results, all approximations show a similarly strong sup¬ 
pression of the signal at large px when using the stated 
fixed reference momentum bin. 

The decorrelation in px observed in Fig. 6 is fairly fast 
and appears incompatible with the experimental obser¬ 
vations. Experimentally only a slow decorrelation can be 
seen in the data when comparing experimental results for 
Vn in p-l-Pb collisions using different methods [63, 66]. 
However a number of caveats are in order with regard 
to this comparison. We emphasize that our results are 
for quarks or more generally on the parton level. While 
hadronization effects will weaken the strong dependence 
on the choice of reference momentum observed on the 
parton level, a quantitative description of the experimen¬ 
tal data in initial state frameworks will also be quite sen¬ 
sitive to the choice of the fragmentation scheme. The role 
of fragmentation in such correlations deserves a more de¬ 
tailed study in the future (see also [67-69]). 

V. THE CGC AND THE COLOR FIELD 
DOMAIN MODEL 

We will now discuss the relation of the azimuthal cor¬ 
relations derived in the CGC framework to those com¬ 
puted recently in the color field domain model introduced 
in [10, 11, 48, 49]. Since the latter qualitatively describes 
some key features of the ridge data in proton-lead colli¬ 


sions at the LHC, it is interesting to compare and con¬ 
trast this model with the CGC based calculations we dis¬ 
cussed thus far. 


A. Electric fields in the Glasma graph 
approximation 

The color field domain model is usually formulated in 
terms of transverse color electric fields and their corre¬ 
lators. To achieve an “apples-to-apples” comparison, we 
will first show how our previous discussion in terms of 
dipole-dipole correlators can be formulated in terms of 
color electric fields. 

The classical color electric fields in the target, as shown 
in Fig. 1, can be expressed in terms of the Wilson line 
correlators as 

Ei{x) = iV{K)d^V^x) . (38) 

Performing a short distance expansion of the dipole op¬ 
erator, one obtains 

D(x,y):^l-^E“(b)E“(b). (39) 

where we denote r = x — y and b = (x -|- y)/2. Within 
the weak field limit of the Glasma graph approximation 
outlined previously, one can evaluate the correlator of 
color electric fields as 

- y), (40) 

where we have used Eq. (21). We emphasize however 
that this equivalence is valid only in the combined limit 
of weak fields and short distances. Using the electric 
field correlator then yields the following expression for 
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the dipole operator 

(X>(x,y)) ~ l + c/-^didj 7(i’)lr=o ■ (41) 

One can similarly use Eq. (39) and express the dipole- 
dipole correlator in the short distance limit as 


(I?(x, y)I?(u, v)) ~ D{^ - y)D{n - v) 


-f 


■ 1 ^ 1 ^ 2'-2 


i?“(bi)E“(bi)E;^(b2)£;f(b2) 


16A^c 

- (£;“(bi)£;“(bi)V£;^(b2)Ef(b2) 


(42) 


B. The color field domain model and non-Gaussian 
correlations 

We focused thus far on conventional models based on 
Gaussian correlations of color fields inside a large nucleus. 
It is now interesting to understand how these relate to the 
color field domain model [10, 11, 48, 49]. In our language, 
the color field domain model is obtained by replacing the 
electric field correlator in Eq. (40) by 


/ \ n 

= -—[S^Hl-A) + 2Aad,] 

xV| 7 (x-y). (46) 


where x = bi -|- ri/2, y = bi — ri/2, u = b 2 -I- r2/2, v = 
b 2 — r2/2. This expression shows that the expectation 
value of the dipole-dipole correlator is sensitive to fluctu¬ 
ations of the color electric fields as characterized by the 
four point correlator. 

In the Glasma graph approximation, the electric field 
is linearly proportional to the charge density and thus 
has Gaussian correlations 

= (43) 

(^El{^)El{y))(^E^,in)El,{^^)) 

+ (£;:(x)E,^(u))(E^'(y)E'(v)) 

+ (E:(x)E;'(v))(if,^u)£;^'(y)). 

The four point correlation function can be expressed in 
terms of the two point function as 

(^El{^)El{y)E^A^)Eliiv)) = (44) 

S‘^^S^‘‘d,d,y{^-y)dkda{u-v) 

- u)djdaiy - v) 
+(5“‘^<5''^a,aa(x - v)a,5fc7(y - u). 

Gombining the above expressions, the dipole-dipole cor¬ 
relator takes the form 

(I?(x, y)I?(u, v)) ~ I4(x - y)7^(u - v) 

+ - ^2)) , (45) 

which agrees precisely with the expansion of the result 
in the double gluon exchange approximation in Eq. (24) 
in the short distance limit jrij ~ |r 2 | <C jbij ~ |b 2 |. 

This simple calculation shows that the physics of fluc¬ 
tuating color electric field domains is implicitly contained 
in the conventional Glasma graph picture. While in the 
short distance (large momentum) limit the dipole-dipole 
correlator is expressed in terms of two and four point cor¬ 
relators of electric fields, there is a one-to-one mapping 
between the statistical properties of these electric fields 
and those of the A’s in the Glasma graph calculation. 


This correlator in the color field domain model depends 
explicitly on the effective degree of polarization A and 
the unit vector a characterizing the direction of the color 
electric field. 

Expectation values of operators within the color field 
domain model are computed by a two step averaging pro¬ 
cedure. In the first step, one performs a Gaussian aver¬ 
age with the modified two point correlation function in 
Eq. (46). The second step consists of an average over all 
possible directions of the chromoelectric fields, a, such 
that 

. (V) 

and 

(a*a%''a'). = --- . 

\ /a g 

Implicit in this two step procedure is a physical as¬ 
sumption about time scales. One assumes that partons 
within a color domain of size ^ 1 /Qs generate a color 
electric field oriented in a particular direction a with a 
likelihood A ranging from 0 — 100 % , which is long lived 
on the time scale of the interaction such that partons in 
the projectile are collimated relative to the direction a of 
this color electric field. A similar picture is implicit in 
the work of [ 12 ] where the net momentum transfer from 
the projectile partons to the target takes the place of the 
color electric field. 

When the life time of the degrees of freedom responsi¬ 
ble for the breaking of rotational symmetry within each 
color domain is much larger than the time scale over 
which one performs the average over the color electric 
fields inside the target in Eq. (46), the orientation of the 
domain appears frozen on that time scale and a separate 
averaging is justifiable. However it is not a priori evident 
that such a separation of time scales exists. In particular, 
it is not clear what would be the intrinsic or dynamical 
scale that separates the time scale over which the color 
electric fields align themselves from the time scale over 
which one performs the average over the different orien¬ 
tations of the field and why such an average should be a 
Gaussian average. 
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Let us now discuss explicitly the calculation of the four 
point correlator of color electric fields in the color field 
domain model. The only correlators that are related to 
physical observables are the ones averaged over all unob¬ 
servable degrees of freedom, including the direction of a. 
Therefore, to understand the correlation structure of the 
model we must compare the two- and four-point func¬ 
tions of the electric field after the full two step average. 

Carrying out the two step averaging procedure for 
two point correlators by averaging Eq. (46) over a us¬ 
ing Eq. (47) one obtains 


((E;:(x)ii;^(y)))_ = V| 7 (x - y). (48) 

The result is independent of the effective polarization A 
and the normalization has been chosen in such a way that 
the two-point function Eq. (48) agrees with the previous 
result in Eq. (40) for a rotationally invariant system in 
the short distance limit^. 

The short distance expansion of the dipole-dipole cor¬ 
relator in Eq. (42) involves the (double) average of the 
four point correlation function which, after the first step, 
can be expressed as 


((E:(x)£;^(y)E,^(u)£;'(v)))_ = 

^(K(x)£;^(y))(E,^u)E;'(v)) 

+ (E:(x)if,"(u))(E;^'(y)E'(v)) 

+ (£;*(x)ii;'(v))(£;,^(u)E^(y))\ . (50) 


Performing the second average with Eq. (48), one obtains 
El{^)El{j)E^,{n)E\{^))) ^ = (51) 


^ab^cd 


gacgbd 


^adgbc 


I \ 

1 - -A'^ (5*^(5'=' -b — -b 


Vr7(x-y)Vr7(u-v) 


(l - ^A^^ -b ^ (<5*^5'=' -b 


X Vt7(x - u)V^7(y - v) 

A^ 


1 - -A^ ] -b — 


X Vt7(x - v)Vr7(y - u). 


^ Note that for a rotationally invariant correlation function 7(r) 
in the short distance limit r ^ 0 we can replace 

djdj7(r)|r=o = —Vr7{r)|r=o. (49) 


Comparing the two point and four point correlation func¬ 
tions of the color electric fields in Eqs. (48) and (51), one 
sees explicitly that the four point correlation function can 
not be expressed in terms of the two point function, i.e. 
the color field domain model is non-Gaussian. Specifi¬ 
cally one finds that the dipole-dipole correlator 


(I?(xi,yi)I?(x2,y2)) D(ri)D(r2) 


■’' 2 )^ - ^Vr7(0)^ 


-b 


ClA^ 

16(Vc2 - 1 ) 


riTo 




(52) 


is a sum of a Gaussian piece (present already in Eq. (45)) 
and non-Gaussian terms proportional to A? induced by 
the two step averaging procedure. The non-Gaussian 
terms are referred to as “disconnected contributions” in 
Ref. [48]. In addition to the small dipole limit r!,rl < 
1/Qs that is assumed in this discussion, one can addi¬ 
tionally take the limit 1/Qs- this case the 

two-gluon exchange approximation becomes exact and 
thus the Glasma graphs and the nonlinear Gaussian are 
equivalent to each other. Even in this limit the .A-terms 
are not suppressed in any way and the color field do¬ 
main model remains different from the other approaches 
considered in this paper. 

It is obvious that the (ri • r 2 )^ terms introduce an ad¬ 
ditional Al-dependent cos 2(f> correlation between the two 
dipoles. Upon Eourier transformation, this modifies the 
angular structure of the correlation between the two pro¬ 
duced particles. What remains unclear at this stage is the 
physical origin of this particular form of non-Gaussian 
correlators as well as the magnitude of the non Gaus- 
sianity characterized by the additional A parameter in 
this model. 


C. Interpretations of the color field domain model 

We noted previously but wish to emphasize again that 
fluctuating domains of color electric field are present also 
in the Glasma graph or non-linear Gaussian approxima¬ 
tion and they are not a new physical feature added by the 
color field domain model. What is different in the color 
field domain model is that the direction of the chromo¬ 
electric field is treated explicitly as a long lived degree of 
freedom. By modifying the correlation function of elec¬ 
tric fields according to Eq. (46) and performing a sepa¬ 
rate average over the orientation of the electric fields the 
statistics of these domains is altered significantly. Most 
importantly, this can lead to sizable non-Gaussian corre¬ 
lations depending on the magnitude of the parameter A 
in this model. Since the single inclusive distribution is 
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not sensitive to A, its value can only be determined from 
correlation measurements. An observable that would be 
particularly sensitive to the presence of intrinsic non- 
Gaussianities would be the four-particle cumulant flow 
coefficient, as discussed in [48]. We will now discuss three 
possible interpretations of the non-Gaussian correlations 
represented by the A-term in the color domain model. 

1. The color electric field is a nonlinear function of 
the color charge density. Thus even if the color 
charges (or, equivalently the A’s) have a Gaussian 
distribution, the color electric fields and thus the 
dipole operators can have non-Gaussian correla¬ 
tions. One possible interpretation of the color field 
domain model is as an effective way to account for 
the non-linear relation between electric fields and 
color charge densities in the target in an otherwise 
linearized calculation. The Glasma graph approxi¬ 
mation of Eq. (38) assumes that the electric field is 
linearly proportional to the color charge and thus 
the electric fields have Gaussian correlations. The 
nonlinear Gaussian approximation, on the other 
hand, has Gaussian correlations for color charges 
but not for the electric fields. In this interpreta¬ 
tion of the color field domain model, the corrections 
encoded in A should be proportional to the differ¬ 
ence between the Glasma graph and the nonlinear 
Gaussian computations in this paper. The total 
anisotropy of the azimuthal two-particle distribu¬ 
tion calculated from the MV model (see e.g. [11]), 
is then the sum of the Glasma graphs and the A- 
term. In this interpretation one has parametri¬ 
cally A^ l/(Nc^ — 1), since correlations are Nc- 
suppressed (relative to the uncorrelated term) in 
both the nonlinear Gaussian and the Glasma graph 
approximations. Our numerical results in Fig. 4 
show that such non-linear corrections can indeed be 
sizeable and should be taken into account as a cor¬ 
rection to the Glasma graph result. However, if this 
is the interpretation it would seem more natural to 
directly use the non-linear Gaussian approximation 
rather than introducing an additional parameter. 
In particular, it is not obvious whether a constant 
A could have a similar momentum dependence as 
the nonlinear Gaussian approximation, given that 
the color field domain model differes from the non¬ 
linear Gaussian even in the small distance limit. 

2. A second possible interpretation of the A-term is 
that it represents non-Gaussian correlations that 
can emerge from JIMWLK rapidity evolution even 
when starting from a Gaussian initial condition. 
This contribution is, in our present calculation, 
represented by the difference between the full 
JIMWLK result and the nonlinear Gaussian. In¬ 
deed we see signs of a ^ 10% deviation between the 
two for pt > Qs- However, this difference is rela¬ 
tively small in practical terms and might not have 
a significant influence on phenomenology. From 


a theory perspective, it is to our knowledge the 
first instance observed in the literature of a mean¬ 
ingful breaking of the Gaussian approximation to 
JIMWLK. We see no obvious reason why the devia¬ 
tion from Gaussianity seems larger here than in the 
observables studied in Ref. [37]. This issue might 
call for additional studies in the future including a 
more systematical check of discretization effects in 
the lattice calculations. 

3. Finally the most intriguing possibility is that the 
A-term represents an intrinsic non-Gaussian cor¬ 
relation that is present in the initial condition for 
JIMWLK evolution and survives substantially af¬ 
ter evolution. This is the interpretation suggested 
in [10, 48]. The possible existence of such non- 
Gaussian correlators was previously suggested in 
[6, 7] and later studied in [70, 71]. While the Gaus¬ 
sian MV model can be justified on quite general 
grounds [72] as arising, due to the central limit the¬ 
orem, from a superposition of a large amount of 
uncorrelated color charges in a heavy nucleus, de¬ 
viations from Gaussian statistics are naturally ex¬ 
pected for a small number of large x degrees of free¬ 
dom. The existence and persistence of such a non- 
Gaussianity at small x would thus be a signal of re¬ 
markably strong long range rapidity correlations in¬ 
side the gluon cascade building up the strong color 
fields at small x. The computations in this paper 
do not address this possibility of an intrinsic non- 
Gaussian four-particle correlation because we have 
been working in the MV model-f JIMWLK evolu¬ 
tion setup where such correlations are absent in the 
initial conditions. 


VI. SUMMARY AND CONCLUSIONS 

We explored a number of calculational schemes to com¬ 
pute two particle correlations of quarks scattering off a 
highly energetic nucleus. All cases correspond to differ¬ 
ent approximations within the dilute-dense limit of the 
color glass condensate framework. The two-particle cor¬ 
relations are quantified in terms of Fourier coefficients in 
an expansion in relative azimuthal angle of the double 
inclusive distribution of scattered quarks. This distribu¬ 
tion is proportional to the dipole-dipole correlator. The 
study of the properties of this correlator in the various 
approximation schemes was the primary objective of this 
work. 

The simplest approximation scheme considered was the 
glasma graph approximation. In this case, the lightlike 
Wilson lines in the dipole-dipole correlator are expanded 
to lowest order, restricting the interaction with the tar¬ 
get to two gluon exchange. One further assumes that the 
gluon correlations in the target are Gaussian correlations. 
This approximation scheme has been used previously in 
the literature to study azimuthally collimated double in- 
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elusive gluon production in p+p and p+Pb collisions. 

Another approximation scheme, of greater complexity, 
is the nonlinear Gaussian approximation. In this case, 
all multigluon exchanges are resummed to all orders to 
obtain a complicated analytical expression for the dipole- 
dipole correlator. This expression is exact as long as there 
are only Gaussian correlations in the target. This is for 
instance the case for the MV model. 

These analytical results are a good benchmark for nu¬ 
merical studies wherein the Wilson lines are computed on 
2-1-1-dimensional lattices for Gaussian distributed sources 
- good agreement is expected and achieved. With the MV 
initial conditions for the Wilson lines at a given rapidity, 
the JIMWLK equations are solved on the lattice to de¬ 
termine the Wilson lines at larger rapidities. These then 
allow one to determine in principle expectation values of 
n-dipole correlators as a function of rapidity. 

We studied the Fourier harmonics V 2 , 'Ca, and V 4 that 
are extracted from azimuthal two particle correlations in 
the various approximation schemes. The Glasma graph 
approximation and the nonlinear Gaussian approxima¬ 
tions differ appreciably for pt ^ 2Qs in the MV model, 
indicating the importance of coherent multiple scatter¬ 
ing effects. Since the Glasma graph approximation is 
at the heart of most comparisons to experimental data, 
this calls for a more detailed study to further quantify 
its theoretical uncertainties. However, we believe that 
in terms of the phenomenological consequences most of 
this difference can be accomodated within the uncertain¬ 
ties in the overall normalization of the Glasma graph 
calculations [1-5]. A significant difference between the 
two approximation schemes is that the symmetry con¬ 
straints inherent in the Glasma graph approximation do 
not produce any odd harmonics in contrast to the non¬ 
linear Gaussian approximation which generates all odd 
harmonics of the azimuthal double inclusive distribution. 
With JIMWLK rapidity evolution, the differences in the 
U2,4 coefficients computed in the Glasma graph and non¬ 
linear Gaussian schemes decrease significantly. Further¬ 
more, we find that the coefficients in the two schemes are 
quite close to those computed by solving JIMWLK nu¬ 
merically without any approximation to the dipole-dipole 
correlator. Since we do not currently have a good inter¬ 
pretation for this better agreement, we believe that it 
may to some extent be accidental. 

We analyzed the dependence of our results on the num¬ 
ber of colors by comparing computations for SU(3) and 
SU(2) gauge helds and found precisely the expected scal¬ 
ing of VnAiPr) with I/(Vc^ — !)• This result confirmed 
that the azimuthal angle dependent correlations are sup¬ 
pressed parametrically by 1 /Nc'^. 

We studied the dependence of the Fourier coefficients 
on the reference transverse momentum (the mo¬ 

mentum of the second scattered quark) in the different 
approximation schemes. These showed clear differences 
for the two reference momenta considered: = px, 

and 0.5 Qs < p^®^ < 3Qg. For all the approximation 
schemes, the choice of equal pr led to a larger signal for 


Pt ^ Qsi with the Glasma graph approximation show¬ 
ing the largest differences. JIMWLK rapidity evolution 
seems to increase the difference between the different p^®^ 
choices. We note however, that the choice of the frag¬ 
mentation scheme can qualitatively influence the com¬ 
parison of model computations of gluon correlations to 
the hadron correlation data. This topic deserves a more 
detailed study in the future. 

Finally we analyzed in detail the relation of the color 
glass condensate based computations to a color do¬ 
main model which captures qualitative features of the 
multiparticle azimuthal correlations observed in proton- 
nucleus collisions. In this framework, the dipole-dipole 
correlator is modified to include an additional term that 
models the polarization of gluon fields in individual do¬ 
mains of color charge within the target. We conclude that 
this term, proportional to the polarization parameter A, 
introduces non-Gaussian correlations amongst the color 
electric fields inside the target nucleus. On the other 
hand, if such non-Gaussianities are not explicitly intro¬ 
duced, the color domain model reduces to the MV model 
in the Glasma graph approximation. 

We also discussed possible origins of non-Gaussian cor¬ 
relations of the color fields of a large nucleus. One pos¬ 
sibility is that these correspond to non-Gaussian correla¬ 
tions induced by JIMWLK rapidity evolution of Gaussian 
correlations at the initial rapidity. However our estimates 
of this effect suggest that such non-Gaussian correlations 
are too small to be relevant phenomenologically. A more 
interesting possibility is that the A polarization term in¬ 
troduced in this model arises from intrinsic four point 
correlations that are significant in the initial condition 
and whose magnitude is preserved with rapidity evolu¬ 
tion. Such correlations would be interesting to study in 
the future. 
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